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In a semi-infinite geometry, a, ID, M-component model of biological evolution realizes micro- 
scopically an inhomogeneous branching process for M — > oo. This implies in particular a size 
distribution exponent r' = 7/4 for avalanches starting at a free end of the evolutionary chain. A 
bulk-like behavior with r' — 3/2 is restored if "conservative" boundary conditions strictly fix to its 
critical, bulk value the average number of species directly involved in an evolutionary avalanche by 
the mutating species located at the chain end. A two-site correlation function exponent tr = 4 
is also calculated exactly in the "dissipative" case, when one of the points is at the border. These 
results, together with accurate numerical determinations of the time recurrence exponent Tf iT3t , 
show also that, no matter whether dissipation is present or not, boundary avalanches have the same 
space and time fractal dimensions as in the bulk, and their distribution exponents obey the basic 
scaling laws holding there. 

PACS numbers: 64.60.Ht,64.60.Ak,05.40.+j,05.70.Jk 



I. INTRODUCTION. 

Nature offers many examples of systems driven by 
some external force towards an out of equilibrium state 
characterized by critical spatiotemporal correlations |Q. 
In this stationary state the accumulated stress is dissi- 
pated by avalanches of activity which occur intermit- 
tently and cover all space and time scales. Models of 
nonequilibrium critical dynamics displaying such features 
have been proposed for several phenomena, ranging from 
earthquakes Q to interface depinning ||, or biological 
evolution in ecosystems 0j. 

Some models of this self-organized criticality (SOC) 
are characterized by extremal dynamics. Among these, 
the model of biological evolution introduced by Bak and 
Sneppen (BS) constitutes an important example 0] . Es- 
pecially for system with extremal dynamics, very few ex- 
act results are available so far ||. Most of our insight 
is based on numerical simulations, scaling arguments ||, 
or mean field solutions 0, related in general to random 
neighbor versions of the models. 

Among the existing mean field approaches, a partic- 
ularly rich and complete one, proposed recently, allows 
to describe in terms of an inhomogeneous branching pro- 
cess (BP) several properties of avalanches, including some 
due to border effects M . In view of its phenomenological 
character, an open interesting problem within such an ap- 
proach remains the identification of specific microscopic 
models realizing the scalings of the inhomogeneous BP in 
some appropriate random neighbor, or similar limit 

A step towards the establishment of an analytical the- 
ory of scaling in extremal dynamics systems was made 



recently by Boettcher and Paczuski who computed 
exactly a correlation function of an M-component ver- 
sion of the BS evolution model in the limit when M ap- 
proaches infinity. This result, combined with numerical 
ones, allowed verification of important general scaling 
laws for self organized critical behavior ||. Such laws 
describe the connection between space and time fractal 
properties of avalanches in the bulk. 

For models like sandpiles, SOC can be established 
thanks to effects of the boundary dissipation, which bal- 
anches the flux of added particles [0. This basic cir- 
cumstance, together with experience with standard crit- 
icality, called recently attention on the boundary scal- 
ing properties of avalanches pLO. Surface scaling in 
sandpiles can be different from the bulk one and can also 
depend on the type of boundary conditions (b.c.) consid- 
ered. In those models the obvious alternative to dissipa- 
tive b.c. are conditions in which part of the border does 
not dissipate grains O]. For evolution models, which do 
not dissipate particles, it is not known whether boundary 
conditions could influence scaling at the border and, in 
case, what should correspond to conservation. These are 
issues we address in the present article. 

Especially in the perspective of obtaining exact results, 
the study of boundary scaling should represent an im- 
portant step towards a deeper and more complete theo- 
retical understanding of SOC. In the present article we 
show that the M — oo model of Ref. [jlO| , if considered 
in the presence of boundaries, provides a microscopic 
realization of the inhomogeneous BP introduced in [g. 
Besides generalizing immediately to this model results 
known for BP, this opens new possibilities of both ana- 
lytical and numerical investigation. In particular, by ex- 
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tending methods used previously for the bulk 10 1, we are 
able to compute exactly the asymptotic two-point cor- 
relator when it involves a point on the boundary in the 
M = oo limit. These results, and an accurate numeri- 
cal analysis of time correlations at the boundary, allow 
to draw a complete scenario of the scalings obeyed by 
boundary and bulk exponents of the system. Boundary 
avalanches have different scaling properties for different 
boundary conditions. However, even in the case of dissi- 
pative border, in which exponents differ from those in the 
bulk, space and time fractal dimensions of the avalanches 
remain unaltered and satisfy the same basic scaling rela- 
tions. 

This paper is organized as follows. In the next sec- 
tion, after introducing the semi-infinite M-component 
model, we discuss its relation with inhomogeneous BP 
in the M — > oo limit and derive a number of analyti- 
cal and numerical results for size and spatial distribution 
properties. The third section is devoted to a discussion 
of exponents related to the time recurrence of activity at 
the border. The last section contains general conclusions. 



II. SEMI-INFINITE BS EVOLUTION CHAIN 
WITH M = oo COMPONENTS. 

The possible relevance of the BS model for evolution, 
as revealed, e.g., by paleonthological records [jl3|, has 
been already discussed in the literature. Here we regard 
this model and its variants as an interesting mathemati- 
cal framework within which SOC dynamics can be stud- 
ied. 

We consider an open chain of species, labelled by an 
index i — 1,2,3... Each species is characterized by M 
independent parameters (traits) Xi a , (a — 1,2, -.,M, 
< Xi a < 1), which somehow quantify the ability of the 
species to survive in connection with M different tasks it 
has to perform in the eco-system. The closer Xi a to 1, 
the more the ability connected to the a-th task, and thus 
the higher the chance that the species avoids mutation. 

The dynamics goes as follows. At every time step the 
smallest Xi a , i.e. the weakest among the traits of all 
species, is identified and replaced by 1. Each one of the 
species which are neighbors along the chain of that, i m im 
with minimum x a , get replaced one of their traits (cho- 
sen at random among the M possible ones) with new 
random numbers extracted independently and with uni- 
form probability in the interval (0, 1). A new minimum 
is then searched and this proceeds so that on long times 
the system self-organizes itself in a stationary state with 
all Xi a uniformly distributed in an interval (A c , 1). 

A A-avalanche is identified with a sequence of muta- 
tions starting at site i m i n with xf . — A, and continuing 
until the current minimum x a remains below A. We call 
s (duration of the avalanche) the total number of minima 
with value below A obtained during the avalanche. 

Rather then considering a translation invariant situa- 
tion as in Ref. jl(|, we think here to a semi-infinite chain, 



with suitable b.c. The probability that in the station- 
ary state a A-avalanche has size s, will thus depend on 
the site j, (j = 1,2,...), where the avalanche started. 
Omitting to explicitate the A-dependence, we indicate 
such probability by Pj(s). The above dynamical rules, 
for M — > oo, lead to write 

Pj(a + 1) = A(l - \)[P j+ i(s) + P^is)} + 

S 

A 2 (*-*') j>l (1) 

s'=0 

Eq. ([!]) is derived on the basis of the same considerations 
made in Ref. |L0|. In first place Eq.(|l|) reflects the fact 
that, with our dynamical rules, the starting active site 
(site j for which xj = A = absolute minimum) can acti- 
vate one of its neighbors with probability A(l— A), or both 
with probability A 2 . In the two cases, of course, one or 
two indipendent avalanches follow, respectively, and the 
global avalanche grows up to a total of s + 1 activated 
sites. The above independence, which allows to write 
Eq.(|l]) in such form, holds in the M — > oo limit, which 
is implicitly assumed here. Indeed, only with M = oo 
the evolution of an avalanche is completely unaffected by 
the fact that a given site has been previously involved 
in the same, or another, avalanche. The effects of this 
kind of conditions can indeed be seen to amount to cor- 
rections of the order i, or higher, in the equations of 
motion. Eq. (Jl|) needs of course to be complemented with 
b.c. These conditions can be written in different ways. 
A first possibility is: 



P 1 (s + l)=XP 2 (s) 



(2) 



This means that, once the boundary site 1 becomes ac- 
tive, it can then activate only site 2 (site does not exist), 
and this occurs with the usual modalities. An alternative 
boundary condition is: 

P 1 (s + 1) = A(l - A)[Px (s) + P 2 (s)} + 

S 

A 2 ^P!( S ')P 2 ( S - S ') (3) 

s'=0 

Eq. (^|) means that for the boundary site, when active, the 
rule of setting x a = 1, valid for j > 1, does not apply. 
On the contrary, j = 1 and j = 2 get now the random 
replacement of one of their traits, as sites j — and j + 1 
in the bulk (see Eq. ([!])). In other words, the role of the 
missing site j — is now played by the site j — 1 itself. 

It is straightforward to recognize that, up to minor 
modifications due to the convention assumed here of re- 
placing by 1 the x associated to the minimum trait, 
Eqs.([y) and (||) have the same structure as those describ- 
ing the inhomogeneous BP in ID of Ref. ||. By intro- 
ducing generating functions Pi(x) — YlV=a Pi( s ) x3 i * 
1 . 2, .., it was found there that Pi(x) ~ 1 + c(l - x) 1 " 7 / 4 , 
for x — > 1~, when A = 1/2. 
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FIG. 1. The log- log plots of the integrated s distribution of 
avalanches starting at site 1, Pi>(s) = P(x)dx with b.c. 
(2) (dots) and (3) (crosses). 

This value of A is such to imply an average number 
2A(1 — A) + 2A 2 = 1 of sites activated by each active site 
in the bulk and coincides with A c [jlOj . The average num- 
ber of sites activated by the border site 1 is instead less 
than unity, according to Eq.(||). All this means that, for 
A = A c = 1/2, Pi ~ s~ T ', with t' = 7/4, when s -> oo. 
This has to be compared with the result ~ s~ 3 / 2 JTc| ] 
holding when the site where the avalanche starts is cho- 
sen in the bulk, and implying a mean field bulk exponent 
r = 3/2§. 

We indicate by N(j,r) the probability that, at A = 
A c = 1/2, an avalanche started at site j never reaches site 
r > j. We are interested in the behavior of N for j = 1 
and large r, which is in turn related to the asymptotic 
radial distribution of avalanches starting at the border of 
the chain. The Markovian nature of avalanche evolution 
allows to write 

N(j, r) = ~ + ~[AT(j + 1, r) + N(j - 1, r)} + 



±N(j + l,r)N(j-l,r), 



for 2 < j < r — 1. The b.c. in Eq.p) imply 



while, since obviously N(r,r) = 0, 



N(r - 1, r) = - + ^N(r - 2, r). 



(4) 



(5) 



(6) 



Let us put now N(j,r) = 1 — f(j,r), so that / repre- 
sents the probability that an avalanche starting in site j 
reaches site r. From the above equations follows: 
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FIG. 2. Log-log plots of P m (r) for b.c. (2) (dots) and b.c. 
(3) (crosses). 

A/(j, r) = - 1, r)f(j + 1, r) 1< j < r - 1 



A/(r-l,r) = i- jj/(r-2,r) 
A/(l,r) = i/(2,r), 



(7) 



where A/(fc,r) is the discrete Laplacian of / at site 
k. Since we are interested in the large r behavior, we 
can pass to a continuum limit, introducing the variable 
z = (J — l)/r. By putting f(J,r) = y(z), we obtain from 
Eq-df 



V"(z) Viz) 



< z < 1, 



(8) 



with boundary conditions y'(l)/r — — 3y(l)/4 + 1/2 and 
y'(0)/r — y(0)/2. Eq. (Q) can be integrated and, after 
some algebra, one finds: 



y(0) ~ 1 

^ 2 1 
- 9^3 



dx 



r(i/3)r(i/6) 



r(i/2) 



20.14. 



(9) 



If the probability that a critical avalanche starting at 
site 1 reaches site r is Pm(r) cx r _r «, we conclude 
t' r = 4, from the fact that f(l,r) oc f r °° Pm (x)dx. The 



above derivation extends the approach of Ref . [|10| , which 
yielded tr = 3 for the bulk Pr(t). 

The presence of a border like that specified by the 
b.c. in Eq.(2), should be expected to make the avalanche 
propagation more difficult compared to the bulk situa- 
tion. Thus, the result t' r = 4 is physically sound com- 
pared to tr = 3. The fact that the radial probability dis- 
tribution function decays with a different exponent when 
the starting point is at the border with b.c. (^), is qual- 
itatively consistent with what we know of the two-point 
correlator at an equilibrium critical point when one of 



the points is fixed at the boundary and not in the bulk 




The exact results above for r' and t' r allow to draw 
a first conclusion on the space fractal dimension, D', of 
avalanches starting at a border with b.c. (2). Assuming 
s oc r D for such an avalanche, leads to D'(1—t') = 1— r' R , 
which follows from Pi(s)ds = Pm(r)dr. Thus, r' = 7/4 
and t' r — 4 imply D' — 4. This dimension coincides with 
the bulk one, D, which satisfies the same kind of relation 
D(l — t) = 1 — tr (!(]]. So, at the boundary, there is no 
distinct space fractal dimension for these avalanches, in 
spite of the different r-exponent. 

We verified numerically the above result for r', by 
simulating the model on open finite chains of length 
N = 200, with M — 100 components. Fig. 1 reports 
our finite size data referring to avalanches starting at the 
border with b.c. given by Eq.(2). The distribution is 
well compatible with the expected r' = 7/4 ( we esti- 
mated t' = 1.78 ± 0.04). Direct simulation allows us 
also to investigate the implications of b.c. (3), which, 
as far as r' is concerned, escape analytical control. For 
these b.c, which keep equal to 1 the avarage number of 
sites activated by the border one, we find numerically 
t' = 1.46 ± 0.04, compatible with r' = r = 3/2 (Fig.l). 
Since the result r' = 7/4 should hold for the inhomoge- 
neous BP as long as 2A(1 — A) + A 2 < 1 ||, we conclude 
that t' = t = 3/2 is peculiar of b.c. (3). Consistently 
one can also show that, with b.c. (3), t' r — tr = 3 ex- 
actly. Thus, also t' r is restored to its bulk value by b.c. 
(3). In the SOC context, similar results were previously 
conjectured, on numerical basis, for the Abelian sandpile 
in two dimensions IOJ]. Indeed, for that model border 
avalanches appeared to possess a toppling distribution 
exponent rather close to the bulk value for conservative 
border, while with boundary dissipation a different r' ap- 
plied Border dissipativity in a BS evolution model 
should then be associated to the fact that the average 
number of sites activated by the extremal site is less than 
the critical bulk value. 

Fig. 2 illustrates a numerical determination of t' r for 
avalanches starting with dissipative b.c. (2). We obtain 
t' r = 4.05 ± 0.05 well consistent with our exact result. 
With b.c. (3) we get t' r — 3.02 ±0.08, compatible in this 
case with the exact t' r =tr = 3 (Fig. 2). 



III. TIME FRACTAL PROPERTIES. 



Avalanches of a BS model posses also time fractal prop- 
erties, revealed, e.g., by the distribution of first return 
times of the activity in a given site (time being mea- 
sured by the number of minima which get replaced during 
the avalanche) . Some general relations among exponents 
connected with time and space fractal behavior in the 
bulk |(| can be easily derived by arguing as follows. 




log(t) 

FIG. 3. Log-log plots of the integrated first return time dis- 
rtibutions in the case of b.c. (2) (dots) and b.c. (3) (crosses). 

If we define as Sfi rs t(t) the probability distribution of 
first return times in a given site, and call n(T) the to- 
tal number of returns in a lapse of time T, we expect 
n{T) oc T d , where d is a time fractal dimension, and 

-^7<x / Sfirst{t)tdt (10) 

Upon putting Sf irst (t) oc t~ T f^ t , we get d = T ftrst - 1. 
Let us then call S a ii(t) oc t~ Tal1 the distribution of times 
for all subsequent returns in a given site. We clearly have 
fx S aU (t)dt oc T d , so that d = l-T aU and Tfi rst +T a u = 2. 
At this point, to link space and time fractal properties is 
sufficient to consider an avalanche as made of a total of 
s activated sites within a (/-dimensional hyperspherical 
region of radius r such that s oc r D . If the avalanche 
has time duration t, we must have s oc r d n(t) = r d n(r z ), 
where z is an exponent connecting space and time (t oc 
r z ). The last relation treats all lattice sites within the 
sphere as equivalent, as far as the return of activity is 
concerned. Now, since in our model s = t by definition, 
z = D also applies. Eventually, one finds: 

d 

T first = 2-— (11) 

This relation was first proposed in Ref. Q for bulk ex- 
ponents. The present derivation at first sight looks ap- 
plicable also to avalanches starting at the border. The 
d = 1, M = oo BS model is an ideal context in which 
to test its validity. In Ref. [jl(^] a numerical estimate of 
t first was obtained which turned out to be well compati- 
ble with the value implied by relation (11) (Tfi rst = 7/4). 
We made a similar determination of Tf irst for avalanches 
starting at both dissipative and nondissipative borders. 
The data are plotted in Fig. 3 , where one can clearly ap- 
preciate that the same values of Tf irst apply in the two 
cases. Indeed, for conservative b.c. (Eq. (3)) we esti- 
mate T first = 1.71 ± 0.03 ~ 7/4. With dissipative b.c. 
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there appears to be a longer transient before the asymp- 
totic time scaling behavior establishes. However, we esti- 
mated r'j irst = 1.71 ±0.04, clearly compatible again with 
7/4. In both cases, of course, the bulk avalanches have a 
distribution with Tfi rs t — 7/4. 

Our results indicate that, like the space fractal dimen- 
sion Z)', also the time fractal dimension dl of boundary 
avalanches is the same as its bulk counterpart, with all 
b.c, and Eq. (11) is always satisfied. 

In Ref. [|| a simulation of the d = 1, M = 1 BS model 
yielded a Tf irst sensibly different from Tfi rs t- If such a 
discrepancy would be confirmed by more systematic and 
asymptotic determinations, one should suspect that the 
identity of D' and D, or even the validity of some scal- 
ing relations like Eq.(ll), are somehow granted here by 
the peculiar, classical character of the M = oo model. 
Anyhow, in such a case, a more complex scaling scenario 
would certainly apply to the model of Ref. || . 



ingredient for the existence of peculiar boundary scaling, 
distinct from the bulk one. Indications that this could 
be a general feature of the SOC state come also from 
previous numerical results for sandpiles |fll"| . 

Our study of the return of activity at the border site 
revealed that dissipativity does not determine a new 
boundary Tj irst exponent, consistent with Eq. (11) and 
with the fact that, like D', for boundary avalanches also 
d! remains unaltered with respect to its bulk value. 

This contrasts, if confirmed by further analysis, awaits 
to be elucidated. 
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IV. CONCLUSIONS. 

The M-component BS model in the limit M — > oo is 
an interesting theoretical laboratory for testing proper- 
ties of the SOC state. With the present work we were 
able to compute analytically in d = 1 the exponents r' 
and t' r referring, respectively, to size distribution and 
space correlation of avalanches starting at a border spec- 
ified by b.c. (2). This extends previous results in Ref. 
fiof , which referred exclusively to bulk properties. In ad- 
dition our formulation allowed to establish a direct link 
between this model and the inhomogeneous branching 
process discussed in Ref. ||. 

The results r' = 7/4 and t' r = 4 show that the space 
fractal dimension D' of border avalanches with b.c. (2) 
remains equal to the bulk one {D 1 = D = 4), in spite 
of the change of these exponents. Complemented with 
numerical tests, these results showed the existence of a 
clear cut distinction between the b.c. in Eq. (2) and 
those expressed by Eq. (3). In analogy with the physics 
of sandpile models, we were led to call b.c. (2) dissipa- 
tive, due to the fact that, in force of them, the boundary 
site, on average, is able to transmit activity to less than 
one site, even if the bulk is critical. This dissipativity is 
responsible for boundary values of the exponents t' and 
t' r different from the bulk ones. On the other hand, when 
b.c. are conservative in the sense speciefied by Eq. (3), 
the existence of a geometrical boundary does not seem 
enough to determine a scaling different from the bulk 
one. 

To border dissipation, which for models like sandpiles 
is a necessary condition for the very establishment of 
the stationary SOC state, could be given here a precise 
meaning also in the context of a model with extremal 
dynamics. In this model dissipation reveals an essential 



[1] P. Bak, C. Thang, and K. Wiesenfeld, Phys. Rev. Lett. 

59, 381 (1987); Phys. Rev. A 38, 364 (1988); D. Dhar, 

Phys. Rev. Lett. 64, 1613 (1990). 
[2] Z. Olami, H. J. S. Feder, and K. Christensen, Phys. Rev. 

Lett. 68, 1244 (1992); K. Ito, Phys. Rev. E 52, 3232 

(1995) . 

[3] K. Sneppen, Phys. Rev. Lett. 69, 3539 (1992). 

[4] P. Bak and K. Sneppen, Phys. Rev. Lett. 71, 4083 (1993); 

[5] See, however, S. Maslov, Phys. Rev. Lett. 77, 1182 

(1996) ; M. Marsili, P. De Los Rios, and S. Maslov, Phys. 
Rev. Lett. 80, 1457 (1998). 

[6] M. Paczuski, S. Maslov, and P. Bak, Phys. Rev. E 53, 
414 (1996). 

[7] H. Flyvbjerg, K. Sneppen and P. Bak, Phys. Rev. Lett. 
71, 4087 (1993); J. de Boer, B. Derrida, H. Flyvbjerg, 
A. D. Jackson, and T. Wettig, Phys. Rev. Lett. 73, 906 
(1994). 

[8] G. Caldarelli, C. Tebaldi, and A. L. Stella, Phys. Rev. 
Lett. 76, 4983 (1996). 

[9] At a purely numerical level, one possible candidate has 
been recently identified in a suitable random neighbor 
model of earthquakes belonging to the family of sand- 
piles 0). See S. Lise and A. L. Stella, Phys. Rev. E 57, 
3633 (1998). 

[10] S. Boettcher and M. Paczuski, Phys. Rev. Lett. 76, 348 

(1996); Phys. Rev. E54, 1082 (1996). 
[11] A. L. Stella, C. Tebaldi, and G. Cardarelli, Phys. Rev. E 

52, 72 (1995). 

[12] E. V. Ivashkevich, D. K. Ktitarev and V. B. Priezzhev, 

J. Phys. A 27, L585 (1994). 
[13] S. J. Gould and N. Eldredge, Nature (London) 366, 223 

(1993). 

[14] K. Binder in Phase Transitions and Critical Phenomena, 
Vol. 8 p. 1, Ed. by C. Domb and J. L. Lebowitz (Aca- 
demic, London, 1983). 



5 



